# ============================================================================
# Factor Loadings, Composite Reliability (CR), and Average Variance Extracted (AVE)
# 计算因子载荷、组合信度和平均方差提取
# ============================================================================

# Load required libraries
library(readxl)      # For reading Excel files
library(lavaan)      # For CFA analysis
library(semTools)    # For CR and AVE calculation

# ============================================================================
# Step 1: Read the Excel file
# ============================================================================

# Read the Excel file
data <- read_excel("C:/Users/wangm/Nutstore/1/AI & Big Data Group/领导底线心智BLM-工作拖延/Supplemental materials/CC/Supervisor Bottom-Line Mentality_Recoded.xlsx")

# Select variables for CFA
cfa_vars <- c("SBL1", "SBL2", "SBL3", "SBL4", 
              "RT1", "RT2", "RT3",
              "GC1", "GC2", "GC3", "GC4", "GC5",
              "JC1", "JC2", "JC3", "JC4", "JC5", "JC6",
              "WP1", "WP2", "WP3", "WP4", "WP5", "WP6", "WP7", "WP8")

# Create subset with only CFA variables
cfa_data <- data[, cfa_vars]

# Convert to numeric
cfa_data <- as.data.frame(lapply(cfa_data, as.numeric))

# Remove rows with missing values
cfa_data <- na.omit(cfa_data)

cat("Sample size:", nrow(cfa_data), "\n\n")

# ============================================================================
# Step 2: Define Five-Factor CFA Model
# ============================================================================

model_5factor <- '
  SBLM =~ SBL1 + SBL2 + SBL3 + SBL4
  RT   =~ RT1 + RT2 + RT3
  WGC  =~ GC1 + GC2 + GC3 + GC4 + GC5
  JC   =~ JC1 + JC2 + JC3 + JC4 + JC5 + JC6
  WP   =~ WP1 + WP2 + WP3 + WP4 + WP5 + WP6 + WP7 + WP8
'

# ============================================================================
# Step 3: Run CFA
# ============================================================================

fit_5factor <- cfa(model_5factor, data = cfa_data, estimator = "MLR")

# ============================================================================
# Step 4: Extract Factor Loadings
# ============================================================================

# Get standardized factor loadings
std_loadings <- standardizedSolution(fit_5factor)
std_loadings <- std_loadings[std_loadings$op == "=~", ]

# Define factor information
factor_info <- list(
  SBLM = list(name = "Supervisor Bottom-Line Mentality", items = c("SBL1", "SBL2", "SBL3", "SBL4")),
  RT   = list(name = "Perceived Red Tape", items = c("RT1", "RT2", "RT3")),
  WGC  = list(name = "Work Goal Clarity", items = c("GC1", "GC2", "GC3", "GC4", "GC5")),
  JC   = list(name = "Job Calling", items = c("JC1", "JC2", "JC3", "JC4", "JC5", "JC6")),
  WP   = list(name = "Work Procrastination", items = c("WP1", "WP2", "WP3", "WP4", "WP5", "WP6", "WP7", "WP8"))
)

# ============================================================================
# Step 5: Calculate CR and AVE
# ============================================================================

# Calculate Composite Reliability (CR)
CR_values <- compRelSEM(fit_5factor, return.total = FALSE)

# Calculate Average Variance Extracted (AVE)
AVE_values <- AVE(fit_5factor)

# ============================================================================
# Step 6: Create Results Table
# ============================================================================

# Create a comprehensive table with all information
results_table <- data.frame()

for (factor in names(factor_info)) {
  factor_name <- factor_info[[factor]]$name
  items <- factor_info[[factor]]$items
  
  for (i in 1:length(items)) {
    item <- items[i]
    
    # Get loading for this item
    loading_row <- std_loadings[std_loadings$lhs == factor & std_loadings$rhs == item, ]
    loading <- if(nrow(loading_row) > 0) loading_row$est.std[1] else NA
    
    # Add row to table
    if (i == 1) {
      # First item: include variable name, CR, and AVE
      results_table <- rbind(
        results_table,
        data.frame(
          Variables = factor_name,
          Items = item,
          Loadings = round(loading, 3),
          CR = round(CR_values[factor], 3),
          AVE = round(AVE_values[factor], 3),
          stringsAsFactors = FALSE
        )
      )
    } else {
      # Subsequent items: only item name and loading
      results_table <- rbind(
        results_table,
        data.frame(
          Variables = "",
          Items = item,
          Loadings = round(loading, 3),
          CR = "",
          AVE = "",
          stringsAsFactors = FALSE
        )
      )
    }
  }
}

# ============================================================================
# Step 7: Display and Save Results
# ============================================================================

cat("=", rep("=", 70), "\n", sep = "")
cat("Factor Loadings, Composite Reliability (CR), and Average Variance Extracted (AVE)\n")
cat("=", rep("=", 70), "\n\n", sep = "")

print(results_table)

cat("\n")
cat("=", rep("=", 70), "\n", sep = "")
cat("Summary by Factor\n")
cat("=", rep("=", 70), "\n\n", sep = "")

# Summary table
summary_table <- data.frame(
  Variables = sapply(factor_info, function(x) x$name),
  CR = round(CR_values[names(factor_info)], 3),
  AVE = round(AVE_values[names(factor_info)], 3)
)

print(summary_table)

# Save results
write.csv(results_table, "Factor_Loadings_CR_AVE.csv", row.names = FALSE)
write.csv(summary_table, "CR_AVE_Summary.csv", row.names = FALSE)

cat("\n")
cat("Results saved to:\n")
cat("  - Factor_Loadings_CR_AVE.csv\n")
cat("  - CR_AVE_Summary.csv\n")
cat("\n")

# ============================================================================
# Step 8: Detailed Factor Loadings (Optional - for verification)
# ============================================================================

cat("=", rep("=", 70), "\n", sep = "")
cat("Detailed Factor Loadings with Standard Errors and Significance\n")
cat("=", rep("=", 70), "\n\n", sep = "")

# Show detailed loadings
detailed_loadings <- std_loadings[, c("lhs", "rhs", "est.std", "se", "z", "pvalue")]
colnames(detailed_loadings) <- c("Factor", "Item", "Loading", "SE", "z", "p-value")
detailed_loadings$Loading <- round(detailed_loadings$Loading, 3)
detailed_loadings$SE <- round(detailed_loadings$SE, 3)
detailed_loadings$z <- round(detailed_loadings$z, 2)
detailed_loadings$`p-value` <- ifelse(detailed_loadings$`p-value` < 0.001, "<0.001", 
                                      round(detailed_loadings$`p-value`, 3))

print(detailed_loadings)
write.csv(detailed_loadings, "Detailed_Factor_Loadings.csv", row.names = FALSE)

cat("\nDetailed loadings saved to: Detailed_Factor_Loadings.csv\n")

